#ifndef _h_model_convex
#define _h_model_convex

#include "geo.h"
using namespace std;

#define forint(i, a, b) for(int i=(int)(a); i<=(int)(b); ++i)

class Convex3D {
	Point *a;
	int n;
	struct Face {
		int i, j, k;
		Face(int i, int j, int k) : i(i), j(j), k(k) {}
		Face() {}
	};
	Face *lst;
	int len;
	Point get_normal(int i, int j, int k) {
		return mult(a[i], a[j], a[k]);
	}
	bool facing(int i, Face f) {
		Point normal=get_normal(f.i, f.j, f.k);
		Point p=a[f.i]-a[i];
		return sgn(dot(normal, p))>0;
	}
	void init() {
		int i=1, j=2;
		while(a[i]==a[j]) ++j;
		int k=1;
		while(sameline(a[i], a[j], a[k])) ++k;
		forint(l, 1, n) {
			forint(t, 0, 1) {
				Face f(i, k, j);
				if(facing(l, f)) {
					len=4;				
					lst[1]=Face(i, j, k);
					lst[2]=Face(j, i, l);
					lst[3]=Face(k, j, l);
					lst[4]=Face(i, k, l);
					return;
				}
				swap(i, j);
			}
		}
		printf("ERRO!\n");for(;;);
	}
	Face *tl;
	unsigned char **cc;
	int tn;
	Point *normal;
	bool same_normal(Point p1, Point p2) {
		double s1=fabs(p1.x)+fabs(p1.y)+fabs(p1.z);
		double s2=fabs(p2.x)+fabs(p2.y)+fabs(p2.z);
		if(sgn(s1)==0) return sgn(s2)==0;
		if(sgn(s2)==0) return sgn(s1)==0;
		if(sgn(s1)!=sgn(s2)) return false;
		p1=p1*s2;
		p2=p2*s1;
		return p1==p2;
	}

	void solve() {
		init();
		for(int i=1; i<=n; ++i) for(int j=1; j<=n; ++j)
			cc[i][j]=0;
		int p=max(n/20,1);
		printf("here\n");
		forint(i, 1, n) {
			if(i%p==0) {printf("computing convex i=%d / %d\n",i,n);fflush(stdout);}
			tn=0;
			for(int j=1; j<=len; )
				if(facing(i, lst[j])) {
					tl[++tn]=lst[j];
					lst[j]=lst[len--];
				} else ++j;
			if(tn==0) continue;
			//printf("#1");
			forint(j, 1, tn) {
				Face f=tl[j];
				++cc[f.i][f.j];
				++cc[f.j][f.k];
				++cc[f.k][f.i];
			}
			//printf("#2");
			forint(j, 1, tn) {
				Face f=tl[j];
				forint(t, 1, 3) {
					if(cc[f.j][f.i]==0)
						lst[++len]=Face(i, f.i, f.j);
					int mm=f.i; f.i=f.j; f.j=f.k; f.k=mm;
				}
			}
			//printf("#3");
			forint(j, 1, tn) {
				Face f=tl[j];
				--cc[f.i][f.j];
				--cc[f.j][f.k];
				--cc[f.k][f.i];
			}
			//printf("#4\n");
		}
	}

public:
	Convex3D(int MAXN) {
		n=MAXN;
		MAXN+=10;
		a=new Point[MAXN];
		lst=new Face[MAXN*5];
		normal=new Point[MAXN*5];
		tl=new Face[MAXN*5];
		cc=new unsigned char*[MAXN];
		for(int i=0; i<MAXN; ++i) cc[i]=new unsigned char[MAXN];
	}
	~Convex3D() {
		delete a;
		delete lst;
		delete normal;
		delete tl;
		for(int i=0; i<n+10; ++i) delete cc[i];
		delete cc;
	}
	void assign(int i, Point p) {
		a[i+1]=p;
	}
	void run() {
		solve();
		for(int i=1; i<=n; ++i)
			normal[i]=Point();
		for(int i=1; i<=len; ++i) {
			Face f=lst[i];
			Point p=get_normal(f.i, f.k, f.j).normalized();
			normal[f.i]+=p;
			normal[f.j]+=p;
			normal[f.k]+=p;
		}
		for(int i=1; i<=n; ++i)
			normal[i].normalize();
	}
	int size() {
		return len;
	}
	Patch get(int i) {
		Face f=lst[i+1];
		return Patch(a[f.i], a[f.j], a[f.k], normal[f.i], normal[f.j], normal[f.k], Color());
		//return Patch(a[f.i], a[f.k], a[f.j], Color());
		//return lst[i+1];
	}
};

#undef forint

#endif